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The Ginzburg-Landau (GL) free energy of crystalline color superconductors is important for understanding 
the nature of the phase transition to the normal quark matter and predicting the preferred crystal structure. So 
far the GL free energy at zero temperature has only been evaluated up to the sixth order in the condensate. To 
give quantitative reliable predictions we need to evaluate the higher-order terms. In this work, we present a new 
derivation of the GL free energy by using the discrete Bloch representation of the fermion field. This derivation 
introduces a simple matrix formalism without any momentum constraint, which may enable us to calculate the 
GL free energy to arbitrary order by using a computer. 
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I. INTRODUCTION 


It is generally believed that the inhomogeneous Larkin- 
Ovchinnikov-Fulde-Ferrell (LOFF) phase appears in a super¬ 
conductor when the pairing between different fermion species 
is under the circumstances of mismatched Fermi surfaces fl- 
i. The mismatched Fermi surfaces are normally induced 
by the Zeeman energy splitting 26p in a magnetic field 0. 
Early studies of the LOFF phase were restricted to ID struc¬ 
tures, which include the Fulde-Ferrell (FF) state with a 
plane-wave order parameter A(z) = Ae 2,qz and the Larkin- 
Ovhinnikov (LO) state 01 with an antipodal-wave order pa¬ 
rameter A(z) = 2A cos(2qz). For ,v-wave pairing at weak cou¬ 
pling, it is known that the FF or LO state exists in a narrow 
window dpi < 5p < 6 p 2 , where the lower critical field dpi = 
0.707Aq and the upper critical field dps — 0.754Ao OHH with 
Aq being the pairing gap at vanishing mismatch. However, 
since the thermodynamic critical field is much lower than dpi 
due to strong orbit effect, it is rather hard to observe the LOFF 
state in ordinary superconductors 0. In recent years, experi¬ 
mental evidences for the LOFF state in some superconducting 
materials have been reported M- 

The studies of color superconductivity in dense quark mat¬ 
ter promoted new interests in the LOFF state Il9l- l2lll . Color 
superconductivity in dense quark matter appears due to_the at¬ 
tractive interactions in certain diquark channels I 22 U 26 I 1 . Un¬ 
der the compact star constraints, different quark flavors ( u , d, 
and s) acquire mismatched Fermi surfaces because of the Beta 
equilibrium and the electric charge neutrality. On the other 
hand, recent experiments on ultracold atomic Fermi gases pro¬ 
vided a controllable way to study the fermion superfluidity 
with population imbalance 1271 12811 . Quark color supercon¬ 
ductors under compact star constraints as well as atomic Fermi 
gases with population imbalance are rather clean systems to 
realize the long-sought LOFF phase. 

In addition to the simple FF and LO states, there exist a 
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variety of crystal structures. The general form of the crystal 
structure of the order parameter can be expressed as ||3tl 

p 

A(r) = Ae 2/9 ”‘ r . (1) 

k=l 


A specific crystal structure corresponds to a multi-wave con¬ 
figuration determined by the P unit vectors n/ c (k = 1,2,..., P). 
Around the tricritical point in the temperature-mismatch 
phase diagram, the LOFF phase can be studied rigorously by 
using the Ginzburg-Laudau (GL) analysis since both the gap 
parameter A and the pair momentum q are vanishingly small 
1311 . It was found that the LO state is preferred near the tricriti¬ 
cal point ll29l - [3lll . However, the real ground state of the LOFF 
phase is still not quite clear due to the limited theoretical ap¬ 
proaches at zero temperature. Various theoretical approaches 
suggested that the LOFF state has a complicated crystal struc¬ 
ture near its phase transition to the normal state llOL 1 1 3l l32l — 
[35ll . A recent self-consistent treatment of the ID modulation 
showed that a solitonic lattice structure is preferred near the 
phase transition to the BCS state (34}]. 

The GL free energy is important for us to understand the 
nature of the phase transition to the normal state and to search 
for the most preferred crystal structure near the phase transi¬ 
tion point. In a pioneer work. Bowers and Rajagopal investi¬ 
gated 23 crystal structures at weak coupling by using the GL 
approach flOll . They evaluated the GL free energy up to the 
sixth order in A, 

Q ° l(A) = Pa A 1 + ~/3A 4 + -yA 6 + 0( A 8 ), (2) 

No 2 H y 


where No is the density of state at the Fermi surface. The 
coefficient a is universal for all crystal structures and is given 
by dot] 


a(dp,q) = -1 + f") - T ln 

2 q \q - dp) 2 


A 2 

A) 


4 (q 2 - dp 2 ) 


■ (3) 


Near the conventional second-order phase transition point 
dp = dp 2 with the optimal pair momentum q = 1.1991 dp, 
we have a - (dp - 6 pf)l 6 p 2 . The GL approach is meaningful 
if the phase transition is of second order or weak first order. 
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The GL coefficients [3 and y for 23 crystal structures were 
evaluated in d. For most of the structures, we find (5 < 0, 
which leads to first-order phase transitions. Among the struc¬ 
tures with y > 0 , the favored one seems to be the body- 
centered cubic (BCC) structure with P = 6 . The order pa¬ 
rameter can be expressed as 

A(r) = 2A [cos(2^) + cos(2 qy) + cos(2gz)]. (4) 


energy to a sufficiently high order in A. Actually, if we can 
evaluate the GL free energy to arbitrary order and determine 
the convergence properties of the GL series, we may even treat 
the strong first-order phase transition within the GL approach. 

The higher-order GL coefficients can be evaluated by using 
the diagrammatic approach used in d- In this approach, to 
evaluate the 2k-th order GL coefficient one needs to sum all 
possible configurations that satisfy the momentum constraint 


Further, it was conjectured that the face-centered cubic (FCC) 
structure with P - 8 is the most preferred structure since 
both [3 and y are negative and their absolute values are the 
largest d- The order parameter can be expressed as 


A(r) = 8 A cos 



cos 




(5) 


For the BCC structure, the GL free energy up to the sixth 
order in A predicts a strong first-order phase transition at 
dp = dp, - 3.6Aq with a large gap parameter A - 0.8Ao 
at dp < dp, d- The prediction of a strong first-order phase 
transition indicates that the GL approach is not valid or the 
higher-order terms in A are important. For the FCC structure, 
the GL free energy up to the sixth order in A actually gives no 
prediction because both fi and y are negative. 

On the other hand, there exist some other approaches that 
do not use the GL approximation. Combescot and Mora em¬ 
ployed Eilenberger’s quasiclassical equation with a Fourier 
expansion for the order parameter 1132113311 . This approach 
predicted that the BCC-normal phase transition is of rather 
weak first order. The upper critical field dju, is only 3.7% 
higher than dpi and A - 0.1 A () at dp < dp, foi l. For the FCC 
structure, it was found that its upper critical field is only 1 . 6 % 
higher than dp 2 and hence it is less favored than BCC. An¬ 
other study employed a solid-state-like approach which cal¬ 
culated the free energy of the BCC structure by directly di¬ 
agonalizing the Hamiltonian matrix in the Bloch space El- 
This approach also predicted that the BCC-normal phase tran¬ 
sition is of weak first order and the upper critical field of BCC 
is slightly higher than dp 2 - 

The contradiction between the predictions from the above 
approaches and from the GL approach indicates that the 
higher-order terms in the GL free energy are rather important 
for a quantitative prediction. Actually, for a first-order phase 
transition, the higher-order expansions are crucial to deter¬ 
mine whether the first-order transition is weak or strong. For 
an intuitive understanding, let us add the eighth-order term in 
the GL free energy for the BCC structure. We have 

= PaA 2 + i/?A 4 + iyA 6 + ^A 8 + <9(A 10 ). (6) 


As a naive example, we assume p > 0. For p — 0, we obtain 
a strong first-order phase transition. However, the first-order 
phase transition becomes weaker and weaker if we increase 
the value of p. For p —> +oo, the phase transition approaches 
second order and dp, —» dpi- On the other hand, if p is small 
or even negative, the higher-order terms such as <9(A 10 ) may 
become important. Therefore, to give more precise predic¬ 
tions within the GL approach we need to evaluate the GL free 


2 k 

2(-D /+ 1 q/ = o (7) 

i= I 

for a set of 2k wave vectors lq 11|2 ■ --<-| 2 a } with q, = qn,. For 
large k, the number of the configurations becomes very large 
for the crystal structures with large number of waves P. More¬ 
over, one needs to introduce 2k Feynman parameters to evalu¬ 
ate the integrals. Therefore, the calculation of the higher-order 
terms is tedious and complicated in this formalism. 

In this work, we introduce a new derivation of the GL free 
energy based on a solid-state physics approach to the crystal 
structures. This derivation provides a simple matrix formal¬ 
ism for the GL coefficients without any momentum constraint. 
We find that the formalism used in |j 10h and our formalism can 
be attributed to two different representations of the fermion 
field: The diagrammatic approach flot] employs the usual con¬ 
tinuous momentum representation, while our derivation uses 
the discrete Bloch representation. Since the matrix operations 
can be easily realized by using a computer, this new formal¬ 
ism may enable us to calculate the GL free energy to arbitrary 
order in A. 

The paper is organized as follows. In Sec. Q]]we briefly 
review the diagrammatic approach to the GL free energy. In 
Sec. Iwe present our new derivation of the GL free energy 
by using the solid-state physics approach. The explicit matrix 
formalism of the GL coefficients for some crystal structures 
are presented in Sec. [TV] We summarize in Sec. [V] The natural 
units h = kv = \ will be used throughout. 


II. CONTINUOUS MOMENTUM REPRESENTATION 

We focus on the general two-flavor pairing at high density, 
low temperature, and weak coupling. In this case, the antipar¬ 
ticle degrees of freedom play no role. Therefore, we can start 
from a general effective Lagrangian for two-flavor pairing at 
high density. The Lagrangian density is given by Etl 

-Ceff = lA + [id, - S( p) +p\fi+ |(^ 1 CT2lA*)(^ T Cr 2 ^), (8) 

where <// = (t// u , <// ( i) T denotes the two-flavor fermion field, 
e(p) is the fermion dispersion with p = —iV, g is a contact 
coupling which represents the attractive interaction, and 0-2 is 
the second Pauli matrix in the flavor space. We shall neglect 
all other internal degrees of freedom, such as color and spin. 
These degrees of freedom contribute a simple degenerate fac¬ 
tor and can be absorbed into the definition of the density of 
state No at the Fermi surface. The fermion chemical poten¬ 
tials are specified by the diagonal matrix p - Amg{p u , pfi) in 
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the flavor space, where 

/r u = A i+5/i, lid = /i - 6ii. (9) 

Here 6/i plays the role of the mismatch between the Fermi 
surfaces. For ultra-relativistic fermion systems such as dense 
quark matter we take e(p) = |p|. The model also works for 
nonrelativistic fermion systems such as cold atomic Fermi 
gases. In this case we take s(p) = p 2 (we set the fermion 
mass M - 1/2 with proper units). Weak coupling requires 
that the pairing gap Aq at dfi - 0 is much smaller than the 
Fermi energy e F - /i. At weak coupling, the physical results 
should be universal if we properly scale the physical quanti¬ 
ties by using the pairing gap Ao and the density of states No at 
the Fermi surface. 

To evaluate the GL free energy, we start from the partition 
function Z of the system. In the imaginary time formalism, it 
is given by 

Z = J [di/i][di// f ]e- s ' s (10) 

with the Euclidean action 



Here r = it is the imaginary time and T is the temperature. 
Once the partition function Z is evaluated, the free energy 
density is given by 


G = --^lnZ. (12) 

Fermionic superconductivity is characterized by a nonzero ex¬ 
pectation value of the difermion fields 

ip(-r, r) = -gifi u (T, r)<Ad(T, r) (13) 

The order parameter A(r) = (ip(T, r)) is static but inhomoge¬ 
neous in the LOFF state. At weak coupling and at low tem¬ 
perature, the order parameter fluctuation becomes negligible 

I 


and we can employ the mean-field approach. With the help of 
the Nambu-Gor’kov (NG) spinor 


T(r,r)= « T ’ r !,. 

\ *Ad ( T ’ r )) 


the mean-field Lagrangian can be expressed as 

|A(r)[ 2 


ZMF = ^(T, T)(-dr - K MF m T , r) - 


where the Hamiltonian operator is given by 
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<u _ I e (P) A(r) 

MF 1 A*(r) -s(p ) +f i-6ny 


(14) 


(15) 


(16) 


For convenience we first consider a finite system in a cu¬ 
bic box defined as x,y,z e [-L/2, L/2] and then set L —» oo. 
Imposing the periodic boundary condition, the fermion mo¬ 
mentum becomes discrete and is given by 


2tt / 

P = y [ tex 


+ me y + nt. 


j, /, m, n • 


(17) 


To convert to the momentum space, we use the Fourier trans¬ 
formation for the fermion fields 


mr) = 


p 0) n 


T' t (r,r) = (18) 


P 0) n 


where V = L 3 is the system volume and a> n — (2n + 1 )nT 
(n e Z) is the fermion Matsubara frequency. For the general 
crystal structure of the order parameter given by (|T|), the mean- 
field action in the momentum space can be evaluated as 


V PA 2 1 \ ^ ^ 

*^MF = TX — ~ ^ ^ (i(O n , p) (iOJ n d(i} n ,Q} n rfip,p' ~ ^(jCOn' ? P X 


(19) 


(i> n ,0J n t p,p 


where p, p' are the discrete momenta given by (fT71) and the 
Hamiltonian matrix F( p p < reads 


/ (f p - 6/u)6 py A(F +) p ,p' \ 

\ A(F_) p , p , (-£ p -5fi)5 p .r }■ 


( 20 ) 


characterizes the crystal structure are given by 


(f ±)p.p- - Z 


p',±qi’ 


k=\ 


( 21 ) 


where q/ c = qn^. 

Here and in the following = e(p) - //. The matrices F ± that The partition function in the mean-field approximation can 
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be evaluated in the momentum representation. We have 


where the matrix F reads 


/' 


Zmf = [d'Hd'rie 


,-Sm 


( 22 ) 


Fpp' - 


0 (F +) p ,p' \ 

(F-) p,p- 0 )■ 


Completing the Gaussian integral, we obtain the free energy Using the derivative expansion, we obtain 





(23) 


£2 — rijsi —A - + 
8 


?EEt*M- 


(j) n 1=1 


where the Trln acts in the momentum space and the NG space. 
The inverse fermion propagator S -1 defined in these spaces is 
given by 


(5 )p p'( ico n ) — ico n 6 p^p' 7 fp p - (24) 


where is the free energy of the normal state. 


T \n 

= _ y Zj Trln 

CO,i 


(Sp 1 ) 

T 



(29) 


(30) 


(31) 


The free energy can be evaluated if we can diagonalize the 
Hamiltonian matrix 77 |lp -. However, this is infeasible because 
the momenta p and p' become continuous in the large vol¬ 
ume limit. We therefore turn to the GL expansion of the free 
energy. The standard field theoretical approach is to use the 
derivative expansion. First, we separate the“BCS” self-energy 
2 A and write 


(S _1 W = (S oV-^aW, (25) 

where the inverse of the free fermion propagator reads 


To obtain the GL free energy, we first complete the trace in 
the NG space. It is easy to show the trace vanishes for odd /. 
After some manipulation, we obtain the GL free energy 


^gl(A) - a^A + ^ 

k=2 


a 2 k A 2 k 

k 


The second-order GL coefficient a -2 is given by 





(32) 


(33) 


<S;V = ( <S *o W (S; ? W ). (20) 

with the matrix elements given by 

( 5 ± 1 )p,p' = +6n + £ p )dp,p'. (27) 

The self-energy term 2 a is given by 

(2 A )p, P ' = ALpp-, (28) 


The higher-order GL coefficients with k > 2 are given by 

a 2k =^J] TT \ (S+F+S - F -> k \- (34) 

D n 

Note that the trace is now taken only in the momentum space. 

Next we complete the trace in the momentum space and 
obtain the integral expressions of the GL coefficients. For the 
second order, we have 


Tr[S+F+S_F_] = 2 „„ 


a=\ b=\ p Pi,p 2 ,P3 
P P 


i(jj„ + 6/u - £ p 


<2=1 b— 1 p p 7 

= ZZi 


1 


1 iOJn + Sfl + £p, 

1 


: V< P 2q„ 


c- >- ^p-p',2q fl . ~ j, 

Ux) n + Ofi — ^ p l(x) n + Of! + £p' 

l l 


3p'-p,-2q 6 


a= 1 p 


ilOn +6/J~^p i(jJ„ + C 5(1+ ^p-2q„ 


(35) 


Therefore, at zero temperature the second-order GL coefficient gg can be expressed as 

_i_i_. 06) 

P g J-oc 2 n J (27 r) 3 iE + 6/j - £ p iE + 6 fi + £ p -2 q 

We notice that a 2 // J is universal for all crystal structures. The above integral suffers from ultraviolet (UV) divergence. In the 
weak coupling limit, the pairing and hence the momentum integral is dominated near the Fermi surface. It is convenient to use 
the regularization scheme that the momentum integral is restricted near the Fermi surface; i.e., -A < |p| — // < A. Using the fact 
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that ()//, q «: A «: /r, we can express o-o as »2 = PNoa( 6 fd, q), where 


a( 6 /d, q) 


gNo 



dp 1 1 

An iE + 6/.1 - g iE + 6 /u +- 2p • q 


(37) 


with Ao = i-i 2 /(2n 2 ) being the density of state at the Fermi surface and p denoting the solid angle. The integrals can now 
be analytically worked out and the cutoff dependence can be removed by using the pairing gap at vanishing mismatch, Ao = 
2Ac 1 /f^ A/ °)_ w e finally obtain the analytical expression given by (0. 

For higher orders ( k > 2), the procedure is similar but becomes tedious. For example, for the fourth order (k = 2), we have 


p p p p 


T,[(S+F+5_F_) 2 ] = £££££ £ 


a= 1 b= 1 c=l d= 1 p P 1 .P 2 .P 7 


() P-P j. C ’P2.P> r 

icO„ + SH ~ f,, ' 1 ' P " 2 " ICO,, + ^ ^ ^ ^ 


P4.P5 s U P(,P7 s 

Op 5 -p 6 ,2q, . x Op 7 _p._2q fl 


P P P P 

a= 1 b= 1 c= 1 


[<,-qf,+q r -q<i 


+ hyU - £p 4 

• 0 Z(tZ 7 T 7 


icOn + 6/o + ^p 6 

1 


ico n + 6/0 - £p icj„ + fyz + ^p_2q„ 


1 


1 


iu)„ + 1 5/1 - ^p-2q a +2q b icOn + 6/0 + ^p~2q„+2q b ~2q c 


(38) 


At weak coupling, the 2k-th order GL coefficient (k > 2) can be generally expressed as 

2k 

& 2 k = No ^ Jik{f\\({ 2 ---(\. 2 k)Sq s fi, q s = ^(-l) ,+1 q,. (39) 

qi.qj.q2* i=i 

Here the summation over each q, means the summation over all P wave vectors q a (a — 1,2,..., P). The quantity .L/fq q 2 ... qa*) 
for general k can be expressed in a compact form. 


£S£*/SDi 


iE + dfj. - £ + 2p • k, iE + 5/u + £ - 2p • 1, ’ 


(40) 


where the momenta k, and 1, are given by (we define qo = 0 
for convenience) 

2i—2 2i-l 

k ' = 2 ( " 1)n+lq ”’ = 2 (-1) ' ,+lq "- (41) 
n= 0 n= 0 

Note that we have set A —> oo since the integral over £ is free 
from UV divergence for k >2. 

The above formalism is completely the same as that ob¬ 
tained by using the diagrammatic approach m. For the non- 
relativistic case, the results can be obtained by replacing the 
pair momentum q with v^q where vf = 2 yfjj is the Fermi 
velocity. The density of state No at the Fermi surface is re¬ 
placed by No = yfjj/(An 2 ). Therefore, at weak coupling, the 
GL free energy is universal for both the relativistic case and 
the nonrelativistic case once we properly express the free en¬ 
ergy in terms of the density of state No at the Fermi surface, 
the pairing gap Ao at vanishing mismatch, and the quantity 
vf q (vf = 1 for the ultra-relativistic case). 

The GL coefficients up to the sixth order ( k < 3) for 23 
crystal structures have been evaluated numerically by Bow¬ 
ers and Rajagopal m. They introduced Feynman param¬ 
eters to evaluate the integral J 2 k(^\^i ■ ■ ■ q 2 *)- The advan¬ 


tage of the above formalism obtained by using the contin¬ 
uous momentum representation is that one can directly ap¬ 
proach the weak coupling limit by using the momentum cut¬ 
off scheme —A < IpI - p < A. However, for higher-order 
coefficients with large k, the calculation becomes complicated 
and tedious. First, for general k, one needs to introduce 2k 
Feynman parameters xi,X 2 and y\,y 2 ,-,yk- At large 
k, the integral over the Feynman parameters becomes com¬ 
plicated. Second, to obtain the GL coefficients, one needs to 
sum over all possible configurations that satisfy the constraint 
£^](-l) ,+ l q,- = 0. At large k, the number of these configura¬ 
tions becomes also large, which makes the calculation tedious. 
To the best of our knowledge, no results of the higher-order 
GL coefficients (k > 4) have been reported so far. 


III. SOLID-STATE PHYSICS APPROACH: DISCRETE 
REPRESENTATION 

In this section, we turn to a discrete representation inspired 
by solid-state physics. For a specific crystal structure given 
by Q. it is periodic in coordinate space. The order parameter 
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can be alternatively expressed as 

A(r) = A/(r), (42) 

where 


which is known as the Bogoliubov-de Gennes (BdG) equa¬ 
tion, is analogous to the Schrodinger equation of a quantum 
particle moving in a periodic potential. The BdG equation for 
the present system can be expressed as 


p 

f(r) = Y J e 2iqikr (43) 

k= 1 

is a periodic function. Here we assume that the function /(r) 
corresponds to a 3D lattice structure. The derivation of the 
GL free energy can be easily generalized to ID and 2D lattice 
structures. For a 3D lattice structure, the unit cell is generated 
by three linearly independent vectors ai, a 2 , and a 3 . We can 
accordingly define the reciprocal space generated by three lin¬ 
early independent vectors bj, b 2 , and b 3 with b, obtained by 
the relation a, • b j = 2n6,j. The periodicity of the order param¬ 
eter means /(r) = /(r + a,). Therefore, the function /(r) can 
be decomposed into a discrete set of Fourier components. We 
have 


/ «(P) -/J-S/j. A(r) 

\ A*(r) -s(p) + ^ - 6 /j 


0,t(r) = E A (f>,i(r). 


(48) 


According to the Bloch theorem, the eigenfunction <p A ( r) takes 
the form of the Bloch function; i.e., 

0,i(r) = e' kl >, lk ( r), (49) 


where the lattice momentum k is restricted in the Brillouin 
zone (BZ) and the function <p A k (r) has the same periodicity 
as the order parameter A(r). We therefore have the similar 
Fourier expansion 


CO 

^k(r) = Yj ^ G e /G r = J] c/>i mn e iG, ™ r . (50) 

G l,m,n =-oo 


/( r ) = ^ fG e ‘ Gr = Y 

G l,m,n=—oo 

oo 

/*( r) = Y& e ~ iGr= E 


Substituting this expansion into the BdG equation, we finally 
obtain a matrix equation in the G-space, 

Y'HgMWg' = E A (k)<p G , (51) 

G' 


where the reciprocal lattice vector G is given by 


where the Hamiltonian matrix < 7 /g,G'( k) is given by 


G = G i mn — /bi + mb 2 + nb 3 , l, m, n e Z. (45) 
The Fourier components /g and f*, are given by 




(( : kl-G - 4/V)4 g,G' A/g_G' 

A/q-_ g (“fk+G - 4p)()G,G' 


■ (52) 


f G = ± fd\f(v)e ' Gr , 

*C Jc 

r a = ±rd 3 rf(r)e iGr , (46) 

Vc Jc 

where V c is the volume of the unit cell and the integration is 
restricted in the cell volume. It is easy to show that 

Y I/gI 2 = P. (47) 

G 

Since the order parameter or pair potential A(r) is periodic, 
the eigenvalue equation for the fermionic excitation spectrum. 


For a given momentum k in the BZ, we can solve the eigen¬ 
values ZGlk) by diagonalizing the above matrix in the discrete 
G-space. This means that the fermionic excitation spectrum 
forms a band structure, in analogy to the energy spectrum of a 
quantum particle moving in a periodic potential. 

Now we turn to the field theory. We consider a finite sys¬ 
tem spanned by three vectors N\ ai, N 2 & 2 , and and as¬ 
sume periodic boundary condition. Then the system con¬ 
tains N 1 N 2 N 3 unit cells and the thermodynamic limit can be 
reached by setting N, —> 00 . In accordance with the Fourier 
expansion (1441 and the matrix equation ( l5ll . we expand the 
fermion field v F(r, r) in terms of the Bloch function rather than 
using the usual momentum representation (IT8l . We write 


T(r,r) = X e ' k r X X qj(lU) "’ k G)e ‘'""' r+ ' G ‘ r ’ 

'V keBZ 0)„ G 

^(r, r) = -1= Y e_!k r Y E k G)e'"" T- ' G r . (53) 

W keBZ w„ G 


One can recover the usual momentum representation ( ITSl by using the fact that k + G can generate all possible momentum p in 
the momentum space. This expansion defines a new representation with two different quantum numbers k and G. We call this 
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Bloch representation. In the Bloch representation, the mean-field action can be evaluated as 

Smf= yj X X X 1?/+((W "’ k ’ G) - h; jn , &v h k , k /K G , G -(k)| k', G'), (54) 


co n ,a) n r k,k'eBZ G,G' 


where the Hamiltonian matrix 7Y(;,G'(k) is given by (l52l >. In 
mathematics, the Bloch representation corresponds to a simi¬ 
larity transformation of the usual momentum representation, 
which makes the Hamiltonian matrix 7Y p p < block diagonal 
with the blocks characterized by the lattice momentum k. 
The functional path integral can be worked out by performing 
integrals over 'V r (ioj n . k. G) and M j (z7jv. k', G') for all possi¬ 
ble values of {/<u„,k, G) and k',G'}. Since the inverse 
fermion propagator is diagonal in the frequency space and the 
k-space, the free energy can be expressed as 


Here the blocks F± defined in the G-space are given by 


(F+) G,G' - /g-G', (F-)g.G’ - Jg’-G- 

Using the derivative expansion, we obtain 


(62) 


p j 

Cl — + —A - + — 

g y 


ZZZy T 'M- 


(63) 


oj„ keBZ Z=1 


where Q\ is the free energy of the normal state. 



T 

V 


Z Z Trln 

keBZ 


(S 'Ig.G' 
T 


(55) 



(64) 


where the Trln acts in the G-space and the NG space. The 
inverse fermion propagator S’ -1 defined in the G-space and 
the NG space is given by 

(S 'Ig.G'O^k, k) = i(i) n SG,G' - 'T^G.G'O 4 ). (56) 


Using the fact that Sq 1 is diagonal and k + G generates all 
continuum momenta p, we recover the usual expression (OTb . 

Completing the trace in the NG space, we find that the trace 
vanishes for odd l. Then the GL free energy can be expressed 
as 


The free energy can be evaluated if we can diagonalize the 
Hamiltonian matrix 'Hg.G' ( k) to obtain all the eigenvalues or 
the band spectrum {^(k)). This is in principle feasible be¬ 
cause the G-space is discrete. Since the matrix has infinite 
dimensions, we should make a truncation —D < Z, m. n < D 
(D e Z + ) to perform the diagonalization. The size of the ma¬ 
trix we need to diagonalize is 2(2 D + l) 3 . To achieve conver¬ 
gence to the limit D — > oo we normally need a large cutoff D, 
which leads to a large computing cost [ l~T5f| . 

Then we turn to the GL expansion based on the expression 
(l55l ) of the free energy. The procedure is the same as we used 
in Sec. II. First, for the inverse fermion propagator S _1 , we 
separate the BCS self-energy and obtain 


Qgl = a 2 A 2 + X ~T A2k - 


k=2 


The second-order GL coefficient is given by 


P T 
a 2 - —I— 

" 8 V 


X X t ^ s + f xS-F-). 

co n keBZ 


The higher-order GL coefficients a 2 k (k > 2) read 

“» = vZZ Tr [< s *' r * s - F ->1' 

CJ„ keBZ 


(65) 


( 66 ) 


(67) 


(S _1 )g,G' = (Og.G'- (Sa)g.G'. < 57 > 

Here the inverse of the free fermion propagator reads 


,c-k _/(^; 1 )g,g' o 

(S° )g,G' 0 (S: l) 


F G.G' = 


(F-) G,G' 0 


(58) 


with the matrix elements given by 

(S ± 1 )G,G' = (ZW;i + 6 p + ^k+G)<5G,G'- (59) 

The self-energy term I A can be expressed as 

(Za)g.G' = AF g ,gs (60) 

where the matrix F is defined as 

0 (F+)g,G' 


(61) 


Note that the trace is now taken only in the G-space. At zero 
temperature and in the large volume limit, we obtain 

P r° dE r d 3 k 


a 2 — —h 

8 


r°° dE r „ 

j (68) 


(2tt) 3 


and 


(* 2 k 


r°° dE r 

J- oo 2/r J BZ 


t/ 3 k 


& 2 k(iE, k). 


) BZ (27r) 3 

for k >2. Here the quantity 3\ 2 k(iE. k) is defined as 


with 


(S±)G,G'(iE. k) = 7 


^G,G' 


iE + 5jj + 4T - c 


(69) 


k(iE, k) = Tr {[S +(/£, k )F + S-(iE, k)F .] k } (70) 


(71) 
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The GL coefficients should be independent of the represen¬ 
tation. Therefore, a 2 / P is still universal for all crystal struc¬ 
ture. We can evaluate it from the FF state. We have 


a 2 = Faf. (72) 

For the FF state we have /(r) = f(z) = e 2lqz , which can be 
regarded as a ID crystal structure with periodicity a = n/q. 
We have the Fourier transformations 

oo oo 

f n e'~\ f(z) = 2 r n e-'~\ (73) 

n=—oo n =—oo 

where the Fourier components /„ and f* are given by 

fn=f:=S n , 1. (74) 

Therefore, the matrices F± read 


( E\ )n,n' — d n — n q Is (F-)tt.ft' ~ d n >— n , 1* 

The GL coefficient o'.( F can be expressed as 

FF i r°° dE r°° k ± dk ± n dk, , , 

of = -+ — nr^(iE,k x , 

g J-oo 2n Jo 2n J 2n 


(75) 


**). (76) 


The trace in J7L can be worked out analytically. We have 


Jl 2 {iE,k ± ,k : ) = Tr[S + F + S-F_] 

Z V ®n,Bl r ^K2,«3 r 

A iE + 5fi - & iE + Sn + u 


- y— 1 --—, 

V + (5// - ;£ + (5ju + 

where is defined as 

£ n (k ± ,k z ) = [k 2 + (k, + 2n<?) 2 ]' -yU. 


(77) 


(78) 


Here v = 1/2 for the ultra-relativistic case and v = 1 for the 
nonrelativistic case. Since k : + 2nq generates all continuous 
momenta in the z direction, the final result for a 2 can be ex¬ 
pressed as (36). 

For higher-order GL coefficients a 2 k (k > 2), the trace 
in the quantity 3\ 2 k{iE, k) becomes tedious. However, be¬ 
cause the G-space is discrete, this procedure can be real¬ 
ized by using a computer with a proper code. Since the 
G-space has infinite dimensions, we first make a truncation 
—D < l,m,n < D to perform the matrix operations. Pre¬ 
cise results can be approached by using a sufficiently large D. 
Some of the GL coefficients may suffer from UV divergence 
and a proper regularization scheme should be implemented. 
For ultra-relativistic dispersion g(p) = |p|, a 2 and or 4 are di¬ 
vergent. For non-relativistic case with e(p) = |p| 2 , only a 2 is 
divergent. The GL free energy can be expressed as 


AA^ = pJA)\y^(AV 

NoSn 2 \ 6 fi) k \ 6 fi) 


(79) 


where the GL coefficients become dimensionless. We have 


a 2 k 


6 n 2k - 2 

No 


a 2 k- 


(80) 


These dimensionless GL coefficients are universal in the weak 
coupling limit. Therefore, we can evaluate them by using 
the non-relativistic dispersion because all the coefficients a 2 k 
(k > 2) are free from UV divergence. On the other hand, in 
this matrix formalism, we need to treat the momenta k and 
G separately. Therefore, the usual momentum cutoff scheme 
-A < £ p < A is not proper to achieve the weak coupling limit. 
The dependence on the chemical potential q becomes explicit 
in the matrix formalism. Weak coupling corresponds to the 
case /r =; £ P » 6 /j. In practice, we can vary the value of ///<)// 
and approach the weak coupling limit. 


IV. FORMALISM FOR SOME CRYSTAL STRUCTURES 

In the final part of this paper, we consider some crystal 
structures of which the order parameters A(r) are real. We 
have f G = f G and hence 

(F+) G,G' = (F-) G,G' = F g , g- = /g-G'- (81) 

In the following, we list the explicit forms of the matrix F, the 
propagators S±, and the GL coefficients a 2 k (k > 2). As we 
pointed out in Sec. [Hi] we shall use the nonrelativistic disper¬ 
sion e(p) = p 2 . In practice, this leads to faster convergence 
and hence is better for numerical calculations. 


A. LO 


The LO state is a superposition of two antipodal plane 
waves (F - 2) with 


At =(0,0,1), n 2 = (0,0,-1). (82) 

The function /(r) can be expressed as 

/(r) = f(z) = 2cos(2qz), (83) 

which forms a ID crystal structure with periodicity a = n/q. 
The Fourier decomposition is given by 


f(z) = J] f n e lniqz , 

n =—oo 

where the Fourier component f„ reads 

fn — d n ,\ dll,— 1- 

The matrix form of the BdG equation is given by 

y < Wn,n'(k)0n'(k) = E A (k)(/> n (k), 

n' 

where the Hamiltonian matrix f H n ,n'{ k) reads 

I (£, - 6 /j) 6 n , n > A \ 

\ Nfn-n' (-& - dju)8 n , n , ) 

with 

£„(kj_, k z ) = k 2 + ( k z + 2 nq ) 1 - yU. 


(84) 

(85) 

( 86 ) 

(87) 

( 88 ) 
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The BZ can be defined as k z e [-q,q]. The GL coefficients 
a 2 k (k > 2) can be expressed as 


&2k 




The matrix elements of S ± and F are given by 


iS ±)n,n' — 


iE + 6/j. + f n ’ 

En,n' — d n y +1 d n y— 1* 

B. Square 


(90) 


The Square state is a superposition of four plane waves (P = 
4) with 

ni= (1,0,0), n 2 = (-1,0,0), 

n 3 = (0,1,0), n 4 = (0,-1,0). (91) 

The function /(r) can be expressed as 

/(r) = f(x, y) = 2 [cos(2^x) + cos(2#y)], (92) 

which forms a 2D crystal structure. The unit cell is generated 
by two linearly independent vectors ai = ae x and a 2 = ae y 
with the lattice spacing a - n/q. The Fourier decomposition 
is given by 


f(x,y)= Yj f^ (lx+my \ 


(93) 


l,m=- 


where the Fourier component //,„ reads 

flm — (^U + h-l)5 m ,0 + 5/,o(^m.l + d m - 1). (94) 

The matrix form of the BdG equation is given by 

T ^[lm],[l'm'](M)<P[l'm'](M) = L,i(k)0[/ m ](k), (95) 


I'm' 


where the Hamiltonian matrix < K[/ m ],[/ , m'](k) reads 


ifIm dfj )()//'(5 mln ' A 

A fl—l',m—m' ( flm d/d)SlJ'dm,n 


(96) 


with 


fimikx, ky, k z ) = ik x + 2Iq) + (k y + 2 mq) + k z - p. (97) 

The BZ can be defined as k x , k y e [-q, q]. The GL coefficients 
a 2 k ik > 2) are given by 

_ r°° dE n dk, r q dky_ r°° dL 
aik J- oo 27r J_ 9 27T J _ 9 In <*> 27r 

Tr [/m] [(S + FS_F)*]. (98) 

The matrix elements of 5 ± and F read 

dl,l'dm,m' 


iS ±)[lm],[I'm'] ~ „ > 

lE + 6/J + fl m 

E[lm],[l'm'] = idl,l '+1 + dl,l'-\)dm,m' 

dij'id mm '+1 + d mm '— i). (99) 


C. BCC 

The BCC state is a superposition of six plane waves (P = 6) 
with 

fij = (1,0,0), n 2 = (-1,0,0), 
n 3 = (0,1,0), n 4 = (0,-1,0), 
n 5 = (0,0,1), n 6 = (0,0,-l), (100) 

The function /(r) can be expressed as 

fir) = 2 [cos(2g.r) + cos(2gy) + cos(2gz)], (101) 

which forms a 3D crystal structure. The unit cell is generated 
by three linearly independent vectors a z = ae x , a 2 = ae v , 
and a 3 = ae : with the lattice spacing a = n/q. The Fourier 
decomposition is given by 


n e 


,2iq(lx+my+nz) 


m= 2 f" 

l,m,n =-oo 

where the Fourier component reads 


( 102 ) 


flmn — if/, I + d/— 1 ) dm,odfi,0 + d/, 0 (d m? j + d m ,-\) &n ,0 

+ d/fidmX) (6 ,,,i + d n - 1) . (103) 

The matrix form of the BdG equation is given by 

^ ] < 7T[/ mH ] ) [/' m ' n '](k)0[/' /K ' tt '](k) — F 4 (k)0[/ mn ](k), (104) 

I'm'n' 

where the Hamiltonian matrix < 7T[; m «],[/'m'n'](k) reads 

if linn d/j)dl,l' d m , m 'S n , n i A f/—l',m—m',n—n' | . I 0 Si 

A fl-l',m-m',n-n' i~flm ~ d/l)6i/6 m , m '6 n y ) 

with 

flmn - ik x + 2 Iq) 1 + iky + 2mq) 2 + (k z + 2nq) 1 - p. (106) 

The BZ can be defined as k x , k y , k z e [-q, q]. The GL coeffi¬ 
cients a 2 k ik > 2) are given by 


(107) 


r°° dE 


C 9 iiA: v 

r c/c 

J oo 2;r v . 

/- 9 2tt _ 

l- 9 2n . 

2tt 


a2k = 

Tr| (.S’+ FS F) k ]. 

The matrix elements of S ± and F read 

&l,V &m,m' &n,n' 


(^±)[/mn],[/'m'nT ~ * ■> 

lE + SfJ, + £Imn 

F[lmn],[I'm'n'] ~ (^/,/'+1 &l,l'—\)&m,m'&n,n' 

^1,1 '1 &m,m'—\)&n,n' 

+ + £«,«'-! )• (108) 
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D. FCC 


The matrix elements of S± and F can be expressed as 


The FCC state is a superposition of eight plane waves (P = 
8) with 

Ai = 4=(U,1). *2 = -A=(—1, —1, —1), 

V3 V3 

«3 = -p(L -1, 1), fl4 = -^=(-l,l,-l), 

V3 V3 

As = -^=(1,-1,-1), As = -U-1,1,1), 

V3 V3 

n 7 = —p(l-1,-1). fig = ^=(-1,-1,1). (109) 
V3 V3 

The function /(r) can be expressed as 


/(, )=8cos i^MlHtii- 


( 110 ) 


which forms a 3D crystal structure. The unit cell is generated 
by three linearly independent vectors a [ = ae x , a 7 = ae y , and 
a 3 = ae z with the lattice spacing a = s/3n/q. The Fourier 
decomposition is given by 


oo 

fix) = 2 f lmn e^ lx+my+nz \ (111) 

l,m,n=—oo 

where the Fourier component reads 

flmn = (d/,1 + d;_i) ((5 m ,i + 6 m - 1 ) (t5„,i + 6 n -\) . (112) 


The matrix form of the BdG equation and the Hamiltonian 
matrix f H[i mn ],[i'm', ,q(k) take the same forms as (104) and (105) 
but with \i, nn given by 




mn 




+ 



2 nq\ 

vf / 


-/i.(113) 


The BZ can be defined as k x ,k y ,k z e [—-^=, -4=]. The GL 
coefficients a 2 k (k > 2) are given by 


r" dE r vi dk x ( 

' v! dk y 

rv! dk z 

J_ 00 2 n 2n J 

-JL 2n „ 

JL 2^ 

V3 

VI 

V3 

Tr[(S + FS_F)*]. 


(114) 


, _ <)/,/- ( \njn’ ^n,n' 

W ±)[lmn],[l'm , n'] — .,. , T. ? 

iE + 6/u + £, m „ 

F [//»«],[/'m'«'] — i&IJ'+l ~^^l,l'—l)i^m,m , + l + l) 

^ i^n,n '+1 - ] )■ (115) 


V. SUMMARY 


In this work, we presented a new derivation of the GL free 
energy of crystalline (color) superconductors and compared it 
with the conventional diagrammatic derivation ifTotl . The di¬ 
agrammatic derivation and our derivation can be attributed to 
the use of two different representations of the fermion field: 
The diagrammatic derivation employs the usual continuous 
momentum representation, while our derivation uses the dis¬ 
crete Bloch representation. In either formalism, we need to 
evaluate the trace of the form Tr ^(5+F + 5 , _F’_) i: J. In the usual 
momentum representation, taking this trace leads to a summa¬ 
tion of all possible configurations that satisfy the momentum 
constraint 0. For large k, the number of these configurations 
becomes also large, which makes the calculation tedious. In 
our formalism, this trace is taken in the discrete Bloch space. 
Therefore, the calculation of the GL coefficients can be com¬ 
puterized: One can generate the matrices F± and S± and per¬ 
form the matrix operations by using a computer with a proper 
code. Once the computerization is realized, we may be able to 
evaluate the GL free energy to arbitrary order in A. With the 
information of the higher-order terms in the GL free energy, 
we can give more reliable predictions for the phase transitions. 
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